#install.packages("pacman")#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)#devtools::install_github("daniel1noble/orchaRd", force = TRUE)library(tidyverse)library(here)library(lme4)library(orchaRd)#library(gptstudio)library(metafor)library(patchwork)library(alluvial)library(ggalluvial)library(easyalluvial)library(ape)library(clubSandwich)library(emmeans)library(MuMIn)library(kableExtra)# making metafor talk to MuMIneval(metafor:::.MuMIn)# install.packages("pak")#pak::pak("MichelNivard/gptstudio")
2 Getting data loaded
Code
#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))#dat_full <- read.csv(here("data/dat_28_06_2023.csv"))dat_full <-read.csv(here("data/dat_19_07_2023.csv"))# add phylogenetic tree - only topologies# TODO? - we could get better tree from birdtree.org# we can do 50 different trees as in # https://academic.oup.com/sysbio/article/68/4/632/5267840tree_top <-read.tree(here("R/birds_MA.tre"))# tree with branch lengthstree <-compute.brlen(tree_top)#plot(tree)# turning it into a correlation matrixcor_tree <-vcv(tree,corr=T)
3 Custom functions
Code
# custom functions#' Title: Contrast name generator#'#' @param name: a vector of character stringscont_gen <-function(name) { combination <-combn(name, 2) name_dat <-t(combination) names <-paste(name_dat[, 1], name_dat[, 2], sep ="-")return(names)}#' @title get_pred1: intercept-less model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred1 <-function (model, mod =" ") { name <-firstup(as.character(stringr::str_replace(row.names(model$beta), mod, ""))) len <-length(name)if (len !=1) { newdata <-matrix(NA, ncol = len, nrow = len)for (i in1:len) { pos <-which(model$X[, i] ==1)[[1]] newdata[, i] <- model$X[pos, ] } pred <- metafor::predict.rma(model, newmods = newdata) }else { pred <- metafor::predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title get_pred2: normal model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred2 <-function (model, mod =" ") { name <-as.factor(str_replace(row.names(model$beta), paste0("relevel", "\\(", mod,", ref = name", "\\)"),"")) len <-length(name)if(len !=1){ newdata <-diag(len) pred <-predict.rma(model, intercept =FALSE, newmods = newdata[,-1]) }else { pred <-predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title mr_results#' @description Function to put results of meta-regression and its contrasts#' @param res1: data frame 1#' @param res1: data frame 2mr_results <-function(res1, res2) { restuls <-tibble(`Fixed effect`=c(as.character(res1$name), cont_gen(res1$name)),Estimate =c(res1$estimate, res2$estimate),`Lower CI [0.025]`=c(res1$lowerCL, res2$lowerCL),`Upper CI [0.975]`=c(res1$upperCL, res2$upperCL),`P value`=c(res1$pval, res2$pval),`Lower PI [0.025]`=c(res1$lowerPR, res2$lowerPR),`Upper PI [0.975]`=c(res1$upperPR, res2$upperPR), )}#' @title all_models#' @description Function to take all possible models and get their results#' @param model: intercept-less model#' @param mod: the name of a moderator all_models <-function(model, mod =" ", type ="homo") {# getting the level names out level_names <-levels(factor(model$data[[mod]])) dat <- model$data#model$data[[mod]] <- factor(model$data[[mod]], ordered = FALSE)# meta-regression: contrasts # helper function to run metafor meta-regression run_rma1 <-function(name) {rma.mv(yi = SMD, V = VCV, mods =~relevel(dat[[mod]], ref = name), random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),test ="t",method ="REML",sparse =TRUE,data = dat) } run_rma2 <-function(name) {rma.mv(yi = SMD, V = VCV, mods =~relevel(dat[[mod]], ref = name), random =list(~1| FocalSpL,~1| RecNo,~gsub("\"", "", mod) | Obs_ID),test ="t",method ="REML",sparse =TRUE,data = model$data) }# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])if (type =="homo"){ model_all <-map(level_names[-length(level_names)], run_rma1) } else { model_all <-map(level_names[-length(level_names)], run_rma2) }# getting estimates from intercept-less models (means for all the groups) res1 <-get_pred1(model, mod = mod)# getting estiamtes from all contrast models res2_pre <-map(model_all, ~get_pred2(.x, mod = mod))# a list of the numbers to take out unnecessary contrasts contra_list <-Map(seq, from=1, to=1:(length(level_names) -1)) res2 <-map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) # creating a table res_tab <-mr_results(res1, res2) %>%kable("html", digits =3) %>%kable_styling("striped", position ="left") %>%scroll_box(width ="100%")# results res_tab}# test #all_models(mod1, mod = "Treat_mod")
4 Preparing data set
Code
# function for calculating varianceVd_func <-function(d, n1, n2, design, r =0.5){# independent designif(design =="among"){ var <- (n1 + n2) / (n1*n2) + d^2/ (2* (n1 + n2 -2)) # variance } else { # dependent design var <-2*(1-r) / n1 + d^2/ (2*(n1 -1)) # variance } var # return variance}# getting Hedges' g - get small size bias corrected effect sizedat_full$SMD <- dat_full$d / (1-3/(4* (dat_full$NTreat + dat_full$Ncontrol) -9))# flipping d dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection# calucating Vddat_full$Vd <-with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))# extra useful function# function for getting mean and sd from median, quartiles and sample size# get_mean_sd <- function(median, q1, q3, n){# sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd# mean <- (median + q1 + q3)/3 # mean# c(mean, sd)# }# observation iddat_full$Obs_ID <-1:nrow(dat_full)dat_full$Phylo <-gsub(" ", "_", dat_full$FocalSpL)# filtering very large variance and also very small sample sizedat_int <- dat_full %>%filter(Vd <10& Ncontrol >2& NTreat >2)#dim(dat_full)#dim(dat_int)# sorting out modality stuff# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) dat_int %>%mutate(Treat_mod =case_when(Treatment =="A"~"A", Treatment =="AV"~"AV", Treatment =="AVG"~"AV", Treatment =="AVM"~"AV", Treatment =="L"~"AVO", Treatment =="O"~"O", Treatment =="OV"~"OV", Treatment =="V"~"V", Treatment =="VG"~"V", Treatment =="VM"~"V", Treatment =="VP"~"V"),# into how manyTreat_No =case_when(Treatment =="A"~1, Treatment =="AV"~2, Treatment =="AVG"~2, Treatment =="AVM"~2, Treatment =="L"~3, Treatment =="O"~1, Treatment =="OV"~2, Treatment =="V"~1, Treatment =="VG"~1, Treatment =="VM"~1, Treatment =="VP"~1),# des it have some add-onsAdd_on =case_when(Treatment =="A"~"No", Treatment =="AV"~"No", Treatment =="AVG"~"Yes", Treatment =="AVM"~"Yes", Treatment =="L"~"No", Treatment =="O"~"No", Treatment =="OV"~"No", Treatment =="V"~"No", Treatment =="VG"~"Yes", Treatment =="VM"~"Yes", Treatment =="VP"~"Yes"), ) -> dat# creating data just for A, V, and AV dat_short <- dat %>%filter(Treat_mod =="A"| Treat_mod =="V"| Treat_mod =="AV")# for add-on, we only need V and AVdat_short_add <- dat %>%filter(Treat_mod =="AV"| Treat_mod =="V")dat <- dat %>%mutate_if(is.character, as.factor)# reordering factors for better visualisation# dat$Treat_mod <- factor(dat$Treat_mod, levels = c("A", "V", "AV", "O", "OV", "AVO"))
5 Exploratory visualization
For Treat_mod (Treatment), we will only visualise A, V, and AV as O $r $, OV, and AVO are much rarer. But for Type (Trait type), we will use all data.
orchard_plot(mod1, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
# result tableall_models(mod1, mod ="Treat_mod")
Fixed effect
Estimate
Lower CI [0.025]
Upper CI [0.975]
P value
Lower PI [0.025]
Upper PI [0.975]
A
0.402
0.211
0.593
0.000
-1.423
2.227
AV
0.472
0.208
0.736
0.000
-1.362
2.307
AVO
0.595
-0.286
1.477
0.185
-1.422
2.613
O
0.258
-0.318
0.833
0.379
-1.646
2.162
OV
0.018
-0.733
0.770
0.962
-1.946
1.983
V
0.432
0.228
0.637
0.000
-1.394
2.259
A-AV
0.070
-0.248
0.389
0.664
-1.772
1.913
A-AVO
0.194
-0.706
1.093
0.673
-1.832
2.219
A-O
-0.144
-0.748
0.460
0.639
-2.057
1.769
A-OV
-0.383
-1.158
0.392
0.332
-2.357
1.590
A-V
0.031
-0.240
0.301
0.825
-1.805
1.866
AV-AVO
0.123
-0.797
1.043
0.793
-1.912
2.158
AV-O
-0.215
-0.847
0.418
0.506
-2.137
1.708
AV-OV
-0.454
-1.250
0.343
0.264
-2.436
1.528
AV-V
-0.040
-0.340
0.260
0.794
-1.880
1.800
AVO-O
-0.338
-1.390
0.714
0.529
-2.436
1.760
AVO-OV
-0.577
-1.711
0.557
0.318
-2.717
1.563
AVO-V
-0.163
-1.067
0.741
0.723
-2.191
1.865
O-OV
-0.239
-1.186
0.707
0.620
-2.286
1.808
O-V
0.175
-0.434
0.784
0.574
-1.740
2.089
OV-V
0.414
-0.364
1.192
0.296
-1.561
2.389
7.2 Treatment with relevant data
Code
VCVs <-vcalc(vi = dat_short$Vd,cluster = dat_short$SubjectID,rho =0.5)mod1a <-rma.mv(yi = SMD,V = VCVs,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Treat_mod, test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at A vs AV and A vs V# TODO - make these hetero model toosummary(mod1a)
mod1b <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~relevel(factor(Treat_mod), ref ="AV"), test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at AV vs V and AV vs A (already done)summary(mod1b)
df AIC BIC AICc logLik LRT pval QE
Full 8 1692.5188 1727.6541 1692.7637 -838.2594 NA
Reduced 6 1707.4108 1733.7623 1707.5532 -847.7054 18.8921 <.0001 NA
Code
orchard_plot(mod1c, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
# result table#res_tab <- all_models(mod1a, mod = "Treat_mod")# res_tab
7.3 Treatment with some additions
Code
# the effect of additions# this is a part of sensitivity analysisVCVs2 <-vcalc(vi = dat_short_add$Vd,cluster = dat_short_add$SubjectID,rho =0.5)mod5 <-rma.mv(yi = SMD, V = VCVs2, random =list(~1|FocalSpL , ~1| RecNo, ~ Treat_mod | Obs_ID), mod =~ Treat_mod*Add_on -1,test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short_add)summary(mod5)
# testing the number of stimulimod4 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Treat_No, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod4)
bubble_plot(mod4,mod ="Treat_No",group ="RecNo",xlab ="The number of simuli",g =TRUE)
7.5 Trait type
Code
# Type of responsesmod2 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2)
# heteroscadasticity model better than the homoscedasticity model# note LifeHistory has lowest variation but this may be expected? # as it is less flexiable (e.g. the number of eggs?)anova(mod2, mod2b)
df AIC BIC AICc logLik LRT pval QE
Full 8 1757.5731 1793.2272 1757.8024 -870.7865 NA
Reduced 6 1785.4590 1812.1996 1785.5923 -886.7295 31.8859 <.0001 NA
Code
orchard_plot(mod2b, mod ="Type",group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)
7.6 Trait categories
Code
# Category of responsesmod3 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Category -1, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod3)
orchard_plot(mod3, mod ="Category",group ="RecNo", xlab ="Standardised mean differnece (SMD)",angle =45,branch.size =3)
7.7 Predactor guild
Code
# Predactor guild# quite heterogeneous# TODO this could be in random effects - think abou thtis a bit latermod6 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ PredGuild -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod6)
orchard_plot(mod11,mod ="ControlType",group ="RecNo",xlab ="Standardised mean differnece (SMD)")
7.13 Sex
Code
# sex# TODO - this could be interesting# what is in males and femalesmod12 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Sex -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod12)
bubble_plot(mod15,mod ="ln_duration",group ="RecNo",xlab ="log(duration in days)",g =TRUE) +geom_point(data = dat,aes(x = ln_duration, y = SMD,color = Type,fill = Type,size =1/sqrt(Vd)), alpha =0.6) +scale_color_discrete() +#+ # how to put the legend for colourguides(color ="legend")
orchard_plot(mod_full,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
orchard_plot(mod_full,mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
int_type <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Type")bubble_plot(int_type, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)
Code
int_trt <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Treat_mod")bubble_plot(int_trt, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)
# displays delta AICc <2candidates_aic2 <-subset(candidates, delta <5) # model averagingmr_averaged_aic2 <-summary(model.avg(candidates, delta <5)) # relative importance of each predictor for all the modelsimportance <-sw(candidates)
dat_fulle <-qdrg(object = mod_fulle, data = dat_short)# marginalized overall mean at vi = 0 and year.c = 0; also weights = "prop" or "cells" average things over proportionally. if not specified, all groups (levels) get the same weights# res_fulle1 <- emmeans(dat_fulle, # specs = ~ sqeffectN,# df = mod_fulle$ddf, # weights = "prop")
Source Code
---title: "**Integration of multimodal cues does not alter mean but reduces variance in avian responses to predators: a systematic review and meta-analysis**"author: "**Kim + Shinichi et al.**"date: "`r Sys.Date()`"format: html: toc: true toc-location: left toc-depth: 3 toc-title: "**Table of Contents**" output-file: "index.html" theme: simplex embed-resources: true code-fold: true code-tools: true number-sections: true #bibliography: ./bib/ref.bib fontsize: "12" max-width: "10" code-overflow: wrapcrossref: fig-title: Figure # (default is "Figure") tbl-title: Table # (default is "Table") title-delim: — # (default is ":") fig-prefix: Fig. # (default is "Figure") tbl-prefix: Tab. # (default is "Table")editor_options: chunk_output_type: consoleexecute: warning: false message: false tidy: true cache: true---## Setting up```{r}#install.packages("pacman")#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)#devtools::install_github("daniel1noble/orchaRd", force = TRUE)library(tidyverse)library(here)library(lme4)library(orchaRd)#library(gptstudio)library(metafor)library(patchwork)library(alluvial)library(ggalluvial)library(easyalluvial)library(ape)library(clubSandwich)library(emmeans)library(MuMIn)library(kableExtra)# making metafor talk to MuMIneval(metafor:::.MuMIn)# install.packages("pak")#pak::pak("MichelNivard/gptstudio")```## Getting data loaded```{r}#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))#dat_full <- read.csv(here("data/dat_28_06_2023.csv"))dat_full <-read.csv(here("data/dat_19_07_2023.csv"))# add phylogenetic tree - only topologies# TODO? - we could get better tree from birdtree.org# we can do 50 different trees as in # https://academic.oup.com/sysbio/article/68/4/632/5267840tree_top <-read.tree(here("R/birds_MA.tre"))# tree with branch lengthstree <-compute.brlen(tree_top)#plot(tree)# turning it into a correlation matrixcor_tree <-vcv(tree,corr=T)```## Custom functions```{r}# custom functions#' Title: Contrast name generator#'#' @param name: a vector of character stringscont_gen <-function(name) { combination <-combn(name, 2) name_dat <-t(combination) names <-paste(name_dat[, 1], name_dat[, 2], sep ="-")return(names)}#' @title get_pred1: intercept-less model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred1 <-function (model, mod =" ") { name <-firstup(as.character(stringr::str_replace(row.names(model$beta), mod, ""))) len <-length(name)if (len !=1) { newdata <-matrix(NA, ncol = len, nrow = len)for (i in1:len) { pos <-which(model$X[, i] ==1)[[1]] newdata[, i] <- model$X[pos, ] } pred <- metafor::predict.rma(model, newmods = newdata) }else { pred <- metafor::predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title get_pred2: normal model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred2 <-function (model, mod =" ") { name <-as.factor(str_replace(row.names(model$beta), paste0("relevel", "\\(", mod,", ref = name", "\\)"),"")) len <-length(name)if(len !=1){ newdata <-diag(len) pred <-predict.rma(model, intercept =FALSE, newmods = newdata[,-1]) }else { pred <-predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title mr_results#' @description Function to put results of meta-regression and its contrasts#' @param res1: data frame 1#' @param res1: data frame 2mr_results <-function(res1, res2) { restuls <-tibble(`Fixed effect`=c(as.character(res1$name), cont_gen(res1$name)),Estimate =c(res1$estimate, res2$estimate),`Lower CI [0.025]`=c(res1$lowerCL, res2$lowerCL),`Upper CI [0.975]`=c(res1$upperCL, res2$upperCL),`P value`=c(res1$pval, res2$pval),`Lower PI [0.025]`=c(res1$lowerPR, res2$lowerPR),`Upper PI [0.975]`=c(res1$upperPR, res2$upperPR), )}#' @title all_models#' @description Function to take all possible models and get their results#' @param model: intercept-less model#' @param mod: the name of a moderator all_models <-function(model, mod =" ", type ="homo") {# getting the level names out level_names <-levels(factor(model$data[[mod]])) dat <- model$data#model$data[[mod]] <- factor(model$data[[mod]], ordered = FALSE)# meta-regression: contrasts # helper function to run metafor meta-regression run_rma1 <-function(name) {rma.mv(yi = SMD, V = VCV, mods =~relevel(dat[[mod]], ref = name), random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),test ="t",method ="REML",sparse =TRUE,data = dat) } run_rma2 <-function(name) {rma.mv(yi = SMD, V = VCV, mods =~relevel(dat[[mod]], ref = name), random =list(~1| FocalSpL,~1| RecNo,~gsub("\"", "", mod) | Obs_ID),test ="t",method ="REML",sparse =TRUE,data = model$data) }# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])if (type =="homo"){ model_all <-map(level_names[-length(level_names)], run_rma1) } else { model_all <-map(level_names[-length(level_names)], run_rma2) }# getting estimates from intercept-less models (means for all the groups) res1 <-get_pred1(model, mod = mod)# getting estiamtes from all contrast models res2_pre <-map(model_all, ~get_pred2(.x, mod = mod))# a list of the numbers to take out unnecessary contrasts contra_list <-Map(seq, from=1, to=1:(length(level_names) -1)) res2 <-map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) # creating a table res_tab <-mr_results(res1, res2) %>%kable("html", digits =3) %>%kable_styling("striped", position ="left") %>%scroll_box(width ="100%")# results res_tab}# test #all_models(mod1, mod = "Treat_mod")```## Preparing data set```{r}# function for calculating varianceVd_func <-function(d, n1, n2, design, r =0.5){# independent designif(design =="among"){ var <- (n1 + n2) / (n1*n2) + d^2/ (2* (n1 + n2 -2)) # variance } else { # dependent design var <-2*(1-r) / n1 + d^2/ (2*(n1 -1)) # variance } var # return variance}# getting Hedges' g - get small size bias corrected effect sizedat_full$SMD <- dat_full$d / (1-3/(4* (dat_full$NTreat + dat_full$Ncontrol) -9))# flipping d dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection# calucating Vddat_full$Vd <-with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))# extra useful function# function for getting mean and sd from median, quartiles and sample size# get_mean_sd <- function(median, q1, q3, n){# sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd# mean <- (median + q1 + q3)/3 # mean# c(mean, sd)# }# observation iddat_full$Obs_ID <-1:nrow(dat_full)dat_full$Phylo <-gsub(" ", "_", dat_full$FocalSpL)# filtering very large variance and also very small sample sizedat_int <- dat_full %>%filter(Vd <10& Ncontrol >2& NTreat >2)#dim(dat_full)#dim(dat_int)# sorting out modality stuff# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) dat_int %>%mutate(Treat_mod =case_when(Treatment =="A"~"A", Treatment =="AV"~"AV", Treatment =="AVG"~"AV", Treatment =="AVM"~"AV", Treatment =="L"~"AVO", Treatment =="O"~"O", Treatment =="OV"~"OV", Treatment =="V"~"V", Treatment =="VG"~"V", Treatment =="VM"~"V", Treatment =="VP"~"V"),# into how manyTreat_No =case_when(Treatment =="A"~1, Treatment =="AV"~2, Treatment =="AVG"~2, Treatment =="AVM"~2, Treatment =="L"~3, Treatment =="O"~1, Treatment =="OV"~2, Treatment =="V"~1, Treatment =="VG"~1, Treatment =="VM"~1, Treatment =="VP"~1),# des it have some add-onsAdd_on =case_when(Treatment =="A"~"No", Treatment =="AV"~"No", Treatment =="AVG"~"Yes", Treatment =="AVM"~"Yes", Treatment =="L"~"No", Treatment =="O"~"No", Treatment =="OV"~"No", Treatment =="V"~"No", Treatment =="VG"~"Yes", Treatment =="VM"~"Yes", Treatment =="VP"~"Yes"), ) -> dat# creating data just for A, V, and AV dat_short <- dat %>%filter(Treat_mod =="A"| Treat_mod =="V"| Treat_mod =="AV")# for add-on, we only need V and AVdat_short_add <- dat %>%filter(Treat_mod =="AV"| Treat_mod =="V")dat <- dat %>%mutate_if(is.character, as.factor)# reordering factors for better visualisation# dat$Treat_mod <- factor(dat$Treat_mod, levels = c("A", "V", "AV", "O", "OV", "AVO"))```## Exploratory visualization For `Treat_mod` (Treatment), we will only visualise `A`, `V`, and `AV` as `O` $r $, `OV`, and `AVO` are much rarer. But for `Type` (Trait type), we will use all data.::: {.panel-tabset}## Treatment vs. Trait types```{r}# reordering dat_shrot$Treat_mod for better visualisationdat_short$Treat_mod <-factor(dat_short$Treat_mod, levels =c("A", "V", "AV"))# Treatment vs Typedat_short %>%group_by(Treat_mod, Type) %>%summarise(n =n()) -> tab#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = Type)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and trait type")```## Treatment vs. Treatment duration```{r}# tuning Treatment duration into a binary variabledat_short %>%mutate(TDration =case_when(duration_days <1~"< 1 day", duration_days >=1~"> 1 day")) -> dat_short# reformatting datadat_short %>%group_by(Treat_mod, TDration) %>%summarise(n =n()) -> tab# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = TDration)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and duration of treatment")```## Treatment vs. Sex```{r}# Treat_mod vs Designdat_short %>%group_by(Treat_mod, Sex) %>%summarise(n =n()) -> tab#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = Sex)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and study setting")```## Treatment vs. Predator type```{r}# Treat_mod vs Designdat_short %>%filter(PredTo !="") %>%group_by(Treat_mod, PredTo) %>%summarise(n =n()) -> tabtab$PredTo <-factor(tab$PredTo, levels =c("A", "N", "B"), labels =c("Adult", "Nestling", "Both"))#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = PredTo)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and study setting")```## Treatment vs. Design```{r}# Treat_mod vs Designdat_short %>%group_by(Treat_mod, Design) %>%summarise(n =n()) -> tab#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = Design)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and study setting")```## Treatment vs. Season```{r}# Treat_mod vs Designdat_short %>%group_by(Treat_mod, Season) %>%summarise(n =n()) -> tab#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = Season)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and study setting")```## Treatment vs. Setting```{r}# Treat_mod vs Designdat_short %>%group_by(Treat_mod, Setting) %>%summarise(n =n()) -> tab#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = Setting)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and study setting")```## Treatment vs Control type```{r}# Treat_mod vs ControlTypedat_short %>%filter(ControlType !="mix") %>%group_by(Treat_mod, ControlType) %>%summarise(n =n()) -> tab#alluvial(tab1[,1:2], freq = tab1$n)# using ggaruvialggplot(tab,aes(y = n,axis1 = Treat_mod,axis2 = ControlType)) +geom_alluvium(aes(fill = Treat_mod)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =6, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Treatment modality and control type")```## Trait type vs. Treatment duration```{r}# Type vs duration_days# turn duration_days into a binary factor: less than 1 and more than 1dat %>%mutate(TDration =case_when(duration_days <1~"< 1 day", duration_days >=1~"> 1 day")) -> datdat %>%group_by(Type, TDration) %>%summarise(n =n()) -> tab2# using ggaruvialggplot(tab2,aes(y = n,axis1 = Type,axis2 = TDration)) +geom_alluvium(aes(fill = Type)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Trait type and duration of treatment")```## Trait type vs. Sex```{r}# Type vs duration_days# turn duration_days into a binary factor: less than 1 and more than 1dat %>%mutate(TDration =case_when(duration_days <1~"< 1 day", duration_days >=1~"> 1 day")) -> datdat %>%group_by(Type, Sex) %>%summarise(n =n()) -> tab# reorderingtab$Sex <-factor(tab$Sex, levels =c("F", "M", "both"), labels =c("Female", "Male", "Both"))# using ggaruvialggplot(tab,aes(y = n,axis1 = Type,axis2 = Sex)) +geom_alluvium(aes(fill = Type)) +geom_stratum(alpha =0.5) +geom_text(stat ="stratum", size =4, aes(label =after_stat(stratum))) +theme(legend.position ="none") +theme(legend.position ="none",axis.text.x =element_blank()) +# remove x-axis labelsylab("Frequency") +xlab("Trait type and duration of treatment")```:::## Meta-analysis### All random effects```{r}##| warning: false# VCV matrixVCV <-vcalc(vi = dat$Vd,cluster = dat$SubjectID,rho =0.5)mod0 <-rma.mv(yi = SMD,V = VCV, random =list(#~1 | Phylo,~1| FocalSpL,~1| RecNo,~1| SubjectID, # incoprated as VCV~1| Obs_ID),# R = list(Phylo = cor_tree),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0)# TODO - think about whether we add this or notrobust(mod0, cluster = dat$SubjectID)round(i2_ml(mod0), 2)orchard_plot(mod0,group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)```### Reduced model```{r}# reduced modelmod0r <-rma.mv(yi = SMD,V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0r)round(i2_ml(mod0r), 2)# comparing two models#anova(mod0, mod0r)```## Meta-regression: uni-moderator (mostly)### Treatmeant with all data```{r}## Treatment - A, V, AV etc #dat$Phylo <- as.character(dat$Phylo)#match(dat$Phylo, colnames(cor_tree))#match(colnames(cor_tree), dat$Phylo)mod1 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Treat_mod -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod1)round(r2_ml(mod1)*100, 2)orchard_plot(mod1, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)# result tableall_models(mod1, mod ="Treat_mod")```### Treatment with relevant data```{r}VCVs <-vcalc(vi = dat_short$Vd,cluster = dat_short$SubjectID,rho =0.5)mod1a <-rma.mv(yi = SMD,V = VCVs,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Treat_mod, test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at A vs AV and A vs V# TODO - make these hetero model toosummary(mod1a)mod1b <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~relevel(factor(Treat_mod), ref ="AV"), test ="t",method ="REML", sparse =TRUE,data = dat_short)# to look at AV vs V and AV vs A (already done)summary(mod1b)orchard_plot(mod1a, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)# modeling heteroscedasticitymod1c <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo, ~ Treat_mod | Obs_ID), mod =~ Treat_mod, test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short)summary(mod1c)# comparision modelsanova(mod1a, mod1c)orchard_plot(mod1c, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)# result table#res_tab <- all_models(mod1a, mod = "Treat_mod")# res_tab```### Treatment with some additions```{r}# the effect of additions# this is a part of sensitivity analysisVCVs2 <-vcalc(vi = dat_short_add$Vd,cluster = dat_short_add$SubjectID,rho =0.5)mod5 <-rma.mv(yi = SMD, V = VCVs2, random =list(~1|FocalSpL , ~1| RecNo, ~ Treat_mod | Obs_ID), mod =~ Treat_mod*Add_on -1,test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short_add)summary(mod5)```### Treatment as an ordinal variable```{r}# testing the number of stimulimod4 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Treat_No, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod4)bubble_plot(mod4,mod ="Treat_No",group ="RecNo",xlab ="The number of simuli",g =TRUE)```### Trait type```{r}# Type of responsesmod2 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2)mod2c <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2c)mod2d <-rma.mv(yi = SMD,V = VCV,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mod =~relevel(Type, ref ="LifeHistory"),test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2d)orchard_plot(mod2,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)round(r2_ml(mod2)*100, 2)# heteroscadasticity modelmod2b <-rma.mv(yi = SMD, V = VCV, mod =~ Type -1, random =list(~1|FocalSpL , ~1| RecNo, ~Type | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod2b)# make other heteromod2e <-rma.mv(yi = SMD,V = VCV,mod =~ Type, random =list(~1|FocalSpL , ~1| RecNo, ~Type | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod2e)mod2f <-rma.mv(yi = SMD,V = VCV,mod =~relevel(Type, ref ="LifeHistory"), random =list(~1|FocalSpL , ~1| RecNo, ~Type | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod2f)# heteroscadasticity model better than the homoscedasticity model# note LifeHistory has lowest variation but this may be expected? # as it is less flexiable (e.g. the number of eggs?)anova(mod2, mod2b)orchard_plot(mod2b, mod ="Type",group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)```### Trait categories```{r}# Category of responsesmod3 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Category -1, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod3)round(r2_ml(mod3)*100, 2)orchard_plot(mod3, mod ="Category",group ="RecNo", xlab ="Standardised mean differnece (SMD)",angle =45,branch.size =3)```### Predactor guild```{r}# Predactor guild# quite heterogeneous# TODO this could be in random effects - think abou thtis a bit latermod6 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ PredGuild -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod6)round(r2_ml(mod6)*100, 2)orchard_plot(mod6, mod ="PredGuild",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Setting```{r}# Settingmod7 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Setting -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod7)round(r2_ml(mod7)*100, 2)orchard_plot(mod7, mod ="Setting",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Season```{r}# Seasonmod8 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Season -1,test ="t",method ="REML",sparse =TRUE,data = dat)mod8b <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Season,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod8b)round(r2_ml(mod8)*100, 2)orchard_plot(mod8,mod ="Season",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Design```{r}# Designmod9 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Design -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod9)mod9b <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Design,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod9b)round(r2_ml(mod9)*100, 2)orchard_plot(mod9,mod ="Design",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Response period```{r}# Response periodmod10 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ ResponsePeriod -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod10)round(r2_ml(mod10)*100, 2)orchard_plot(mod10,mod ="ResponsePeriod",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Control type```{r}# control typemod11 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ ControlType -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod11)round(r2_ml(mod11)*100, 2)orchard_plot(mod11,mod ="ControlType",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Sex```{r}# sex# TODO - this could be interesting# what is in males and femalesmod12 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Sex -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod12)round(r2_ml(mod12)*100, 2)orchard_plot(mod12,mod ="Sex",group ="RecNo",xlab ="Standardised mean differnece (SMD)")# shoter data for just males and females# hetero but no sex effect heredat_sex <- dat %>%filter(Sex !="both")VCV3 <-vcalc(vi = dat_sex$Vd,cluster = dat_sex$SubjectID,rho =0.5)mod12a <-rma.mv(yi = SMD, V = VCV3, mod =~ Sex, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), #struct = "DIAG",test ="t",method ="REML", sparse =TRUE,data = dat_sex)mod12b <-rma.mv(yi = SMD, V = VCV3, mod =~ Sex, random =list(~1|FocalSpL , ~1| RecNo, ~Sex | Obs_ID), struct ="DIAG",test ="t",method ="REML", sparse =TRUE,data = dat_sex)summary(mod12b)anova(mod12a, mod12b)orchard_plot(mod12b,mod ="Sex",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Age```{r}# agemod13 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Age -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod13)round(r2_ml(mod13)*100, 2)orchard_plot(mod13,mod ="Age",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Predactor type```{r}# type of predatordat$PredTo[dat$PredTo ==""] <-NAmod14 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ PredTo -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod14)round(r2_ml(mod14)*100, 2)orchard_plot(mod14,mod ="PredTo",group ="RecNo",xlab ="Standardised mean differnece (SMD)")```### Treatment duration```{r}# treatment durationdat$ln_duration <-log(dat$duration_days)mod15 <-rma.mv(yi = SMD,V = VCV,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mods =~ ln_duration,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod15)mod16 <-rma.mv(yi = SMD,V = VCV,random =list(~1|FocalSpL,~1| RecNo,~1| Obs_ID),mods =~ ln_duration*Type,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod16)round(r2_ml(mod15)*100, 2)bubble_plot(mod15,mod ="ln_duration",group ="RecNo",xlab ="log(duration in days)",g =TRUE) +geom_point(data = dat,aes(x = ln_duration, y = SMD,color = Type,fill = Type,size =1/sqrt(Vd)), alpha =0.6) +scale_color_discrete() +#+ # how to put the legend for colourguides(color ="legend")#p + geom_point(aes(colour = Type))#scale_colour_manual(values = c("red", "blue", "green"))```## Meta-regression: multi-moderator (mostly)### full models```{r}######################## Mulit-variable models#######################dat_short$sln_duration <-scale(log(dat_short$duration_days))mod_full <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~#Design + sln_duration*Type + sln_duration*Treat_mod + Sex,test ="t",method ="REML",sparse =TRUE,data = dat_short)summary(mod_full)round(r2_ml(mod_full)*100, 2)orchard_plot(mod_full,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)orchard_plot(mod_full,mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)int_type <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Type")bubble_plot(int_type, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)int_trt <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Treat_mod")bubble_plot(int_trt, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)``````{r}# mulit-model selectioncandidates <-dredge(mod_full, trace =2)# displays delta AICc <2candidates_aic2 <-subset(candidates, delta <5) # model averagingmr_averaged_aic2 <-summary(model.avg(candidates, delta <5)) # relative importance of each predictor for all the modelsimportance <-sw(candidates)```## Publication bias### Funnel plot: uni-moderator```{r}```### Funnel plot: multi-moderator```{r}funnel(mod0r, yaxis="seinv",type ="rstudent")``````{r}funnel(mod_full, yaxis="seinv",type ="rstudent")```### Egger regression: uni-moderator```{r}# Eggerdat$effectN <- (dat$Ncontrol * dat$NTreat) / (dat$Ncontrol + dat$NTreat)dat$sqeffectN <-sqrt(dat$effectN)mod0e <-rma.mv(yi = SMD,V = VCV,mods =~ sqeffectN,random =list(#~1 | Phylo,~1| FocalSpL,~1| RecNo,#~1 | SubjectID, # incoprated as VCV~1| Obs_ID),#R = list(Phylo = cor_tree),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0e)bubble_plot(mod0e,mod ="sqeffectN",group ="RecNo",xlab ="Effective N",g =TRUE)```### Decline effect: uni-moderator```{r}# decline effectmod0d <-rma.mv(yi = SMD,V = VCV,mods =~ Year,random =list(#~1 | Phylo,~1| FocalSpL,~1| RecNo,#~1 | SubjectID, # incoprated as VCV~1| Obs_ID),#R = list(Phylo = cor_tree),test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod0d)bubble_plot(mod0d,mod ="Year",group ="RecNo",xlab ="Publication year",g =TRUE)```### All together```{r}# full modeldat_short$effectN <- (dat_short$Ncontrol * dat_short$NTreat) / (dat_short$Ncontrol + dat_short$NTreat)dat_short$sqeffectN <-sqrt(dat_short$effectN)mod_fulle <-rma.mv(yi = SMD, V = VCVs, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ sqeffectN + Year + sln_duration*Type + sln_duration*Treat_mod + Sex,test ="t",method ="REML",sparse =TRUE,data = dat_short)summary(mod_fulle)dat_fulle <-qdrg(object = mod_fulle, data = dat_short)# marginalized overall mean at vi = 0 and year.c = 0; also weights = "prop" or "cells" average things over proportionally. if not specified, all groups (levels) get the same weights# res_fulle1 <- emmeans(dat_fulle, # specs = ~ sqeffectN,# df = mod_fulle$ddf, # weights = "prop")```
